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It is demonstrated that the canonical distribution for a subsystem of a closed system follows directly from the 
solution of the time-reversible Newtonian equation of motion in which the total energy is strictly conserved. It 
is shown that this conclusion holds for both integrable or nonintegrable systems even though the whole system 
may contain as little as a few thousand particles. In other words, we demonstrate that the canonical distribution 
holds for subsystems of experimentally relevant sizes and observation times. 

PACS numbers: 05.20.-y, 05.45.Pq, 75.10.Pq 



Boltzmann's postulate of equal a-priori probability is the 
cornerstone of statistical mechanics. This postulate elimi- 
nates the difficulties of deriving the statistical properties of a 
many-body system from its dynamical evolution by introduc- 
ing a probabilistic description in terms of the (micro)canonical 
ensemble. In classical mechanics, the equivalence of these 
two fundamentally different levels of description remains elu- 
sive HH. 

Although considerable progress has been made through the 
development of ergodic theory, more than hundred years after 
its conception a direct demonstration that the equal a-priori 
principle or, equivalently, the canonical distribution follows 
directly from the dynamical behavior of a finite number of 
particles is still missing QIQ]- 

The discovery of "deterministic chaos" has made a huge 
impact on old discussions of the relations between classi- 
cal and statistical mechanics OHa]. The irreversibility in the 
macroworld is often related with unstable motion in generic 
mechanical systems characterized by positive Lyapunov ex- 
ponents and a non-zero Kolmogorov entropy. Integrable sys- 
tems are considered from this point of view as exceptional 
(not to say pathological) and irrelevant for the problem of 



justification of statistical physics. The famous Fermi-Pasta- 
Ulam paradox H-QxJ and its interpretation in terms of close- 
ness of their model to the completely integrable KdV model 
(see Ref. 0]) has emphasized (probably, overemphasized) this 
point. Indeed, there is no tendency to equilibration in a sys- 
tem of noninteracting entities (particles of an ideal gas, nor- 
mal modes in harmonic oscillator systems, solitons in KdV 
systems, etc.). However, if we discuss an isolated Hamilto- 
nian system there is no way to obtain the statistical mechanical 
behaviour (in particular irreversibility), not even for a system 
with chaotic motion. An alternative is to consider an open sys- 
tem. For the open system, the role of integrability should be 
discussed in a different context. The integrability of the iso- 
lated Hamiltonian system implies that its Liouville operator L 
is diagonal in the representation of angle-action variables Jj§]. 
If we choose a subsystem S of the isolated integrable system 
which is also integrable, then its Liouville operator L$ is di- 
agonal in the angle-action variables of the subsystem. How- 
ever, these variables can be different from those of the isolated 
Hamiltonian system. L and L$ do not commute and cannot be 
diagonalized simultaneously. It is not clear a-priori whether 
in this situation the subsystem S can equilibrate or not, and 



what the conditions for this equilibration are. Surprisingly, 
this very natural issue is not clarified yet and thus requires 
additional studies. 

Here we present clear and unambiguous evidence that the 
canonical distribution for S follows directly from the solution 
of the time-reversible Newtonian equation of motion. Fur- 
thermore, it is shown that the energy of the subsystem and 
of the environment equilibrate on a relatively short, micro- 
scopic time scale, even if the whole system contains only of 
the order of a thousand degrees of freedom. The key feature 
of our demonstration is that we follow the time evolution of 
a closed (isolated) system with a fixed energy and consider 
only a subsystem S of the closed system. The time evolu- 
tion of the entire system is obtained by solving the Newtonian 
equations of motion for typical initial states. We do not per- 
form ensemble averaging nor do we invoke arguments based 
on (non)ergodicity yj, |2j, 1 1 1 1 - 4 1 3fl or on the thermodynamic 



limit QJJ, |2y. Nevertheless, we observe that S is governed by 
the canonical distribution. 

Our demonstration is not a mathematical proof but is based 
on the exact numerical solution of what is perhaps the sim- 
plest of all interacting many-body systems: a one-dimensional 
harmonic oscillator model of a solid. By solving the Newto- 
nian equation of motion of the whole, integrable system and 
analyzing only a single trajectory of the subsystem in phase 
space, we show that the number of times that the subsystem is 
observed to posses a certain energy is distributed according to 
canonical ensemble theory, implyin g fo r example that within 
the subsystem equipartition holds Il44l7ll . Repeating the 
analysis for a one-dimensional model of classical magnetic 
moments which is known to exhibit chaotic behavior 11811 . it 
is found that this conclusion does not depend on whether or 
not the motion is chaotic. For both models, the distributions 
extracted from the single-trajectory Newtonian dynamics are 
found to be in excellent agreement with the corresponding re- 
sults of microcanonical Monte Carlo simulations. 

As a first example, we consider the most basic classical 
model for the vibration in a solid, namely a set of particles of 
mass m, arranged on a ring and connected to their two neigh- 
bors by harmonic springs, see Fig.Q] The Hamiltonian of the 
system reads 
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where Xj and pi are the displacement and momentum of the 
rth oscillator, mQ. 2 is the spring constant and N is the total 
number of particles. For convenience, m and Q. are set to 
1 and all quantities such as momenta and displacements are 
taken to be dimensionless. The positions are constrained to lie 
on a circle. Changing to normal-mode coordinates {P^,^}, 
the Hamiltonian reads H = (1/2) Zk=o { P k +(0 k X k) where 
(Ok = 2\smnk/N\. The motion of each set of coordinates 
{Pk,Xk} is described by a single sinusoidal oscillation and is 
decoupled from the motion of all other sets. Clearly, this clas- 
sical system is integrable which allows us to compute numer- 
ically the coordinates and momenta of the oscillators without 




FIG. 1 . Picture of the harmonic oscillator (large) and magnetic mo- 
ment (small) models, subject to periodic boundary conditions. Parti- 
cles are connected by harmonic springs or carry a magnetic moment 
that interacts with its nearest neighbors. The blue (red) colored par- 
ticles belong to the subsystem (environment). 



introducing systematic or cumulative errors. To this end, we 
transform the set of values of {x\ , . . . ,xn,pi ,Pn} at time 
t = to their corresponding normal-mode values, use the sim- 
ple sinusoidal dependence of the latter to obtain their values at 
any point t in time, and use the inverse transformation to find 
the values of {xi, . . . ,xn, pi, - ■ ■ ,pn} at time t. This whole 
procedure is numerically exact, up to machine precision. 

As a second example, we consider a ring of classical 
magnetic moments, their total energy given by the Hamilto- 
nian 111 



ff = -/£SrS i+ i, 

i=i 



(2) 



where S, is a 3-dimensional unit vector, representing the mag- 
netic moment of a particle at lattice site i, J defines the energy 
scale which we set equal to 1 in our numerical work and N is 
the total number of moments. The equation of motion of these 
moments reads 



-7SiX(Si_i + S /+ i), 



(3) 



which obviously is nonlinear. Nevertheless, Eq. (O admits 
a harmonic-wave solution S,(f) = (acos 0,- + bsin0,)cos + 
csin0, where 0, = ip — (Ot, (0 = 2(1 — cos p) sin0, and p 
are real constants, and (a.b.c) form a right-handed set of or- 
thogonal unit vectors O20l 12111 . More generally, Eq. (O has 
simple analytical solutions for N ~2 and N = 3 11811 . The 
motion of N = 4 magnetic moments arranged on a ring is reg- 



ular 111 8fl . For N > 4, the system exhibits chaotic motion 11811 . 
except for special initial conditions such as the spin-wave and 
soli ton solutions l22ll . We integrate the nonlinear equations of 
motion using a fourth-order Suzuki-Trotter product-formula 
method which conserves (1) the volume of the phase space, 
(2) the length of each magnetic moment and (3) the total en- 
ergy 12311 . Due to the chaotic character of the motion of the 
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FIG. 2. Left: Poincare-type map of the differences AX = (xq — x\ ) /2 
and AP = (po — p\)j1\f2 for a ring of N = 512 oscillators. This 
system is integrable. Right: Poincare-type map of the differences 
AX' = S \ - S\ and AP' = S\ - S§ for a ring of N = 512 magnetic 
moments. This system is chaotic flUl . Each graph contains 5000 
points. 



magnetic moments, their trajectories are unstable with respect 
to rounding and time-integration errors. Nevertheless, the nu- 
merical method used guarantees that the motion of the mag- 
netic moments strictly conserves the energy, as required for a 
microcanonical ensemble simulation. 

In our numerical work, the initial values of the coordinates 
or magnetic moments are generated by standard Monte Carlo 
methods 12411 . see the Supplementary Information. 

In the normal-mode representation, each of the oscillators 
traces out an ellipse in the (X^P^) plane. However, in the 
original coordinates, this simplicity is lost as illustrated by 
the Poincare map shown in Fig.|2|left). Clearly, it is difficult 
to detect some regularity in this map. Moreover, this map 
is very similar to Fig. fright), the Poincare-type map of a 
one-dimensional classical model of magnetic moments. This 
similarity exists in spite of the fact that the oscillator system 
is not chaotic. 

The whole system is divided into two parts, a subsystem 
with Ns particles and an environment with Ne =N — Ns par- 
ticles. The Hamiltonian is written as H = H$ + He + Hse, 
where H$ (He) denotes the energy of the subsystem (environ- 
ment) and Hse denotes the energy due to the interaction of the 
subsystem with the environment. Thus, in the case of particles 
arranged on a ring, H$ and He are open chains of particles, 
and Hse only contains two terms of the form (x, — x i+ i) 2 or 
S; • Si+i for the oscillator and magnetic system, respectively. 

We first study the dynamic evolution of the subsystem when 
it is brought in contact with the environment. Initially, using 
one of the procedures described in the Supplementary Infor- 
mation, the subsystem and environment are prepared such that 
they have a different energy. As the whole system evolves in 
time, strictly conserving the total energy, we monitor the en- 
ergy of the subsystem as a function of time. Some represen- 
tative results are shown in Fig. [3] for both the harmonic os- 
cillator and magnetic moment model. From Fig.[3j it is clear 
that for both models, the energy of the subsystem E s rapidly 
approaches the average energy of the whole system. 

Having established that the Newtonian equation of motion 
drives the subsystem and environment to a common equilib- 
rium state, the next step is to study the distribution of the sub- 




FIG. 3. Time evolution of the energy per particle Es /Ns, as obtained 
from a chain of Ns = 80 particles embedded in a ring of N = 5 1 2 par- 
ticles. Initially, the subsystem of Ns = 80 particles is put in its ground 
state and the configuration of the environment of Ne = N — Ns = 432 
particles is drawn randomly from its canonical distribution at temper- 
ature T = 1. Red line: energy per particle of the subsystem; horizon- 
tal black line: energy per particle of the whole system. Main figure: 
the energy of the oscillators in the subsystem quickly relaxes to the 
average energy of the whole system. The inset shows the correspond- 
ing data for the spin system. 



system energy. After the relaxation to the equilibrium state, 
we monitor the energy Es of the subsystem at regular time 
intervals and construct a (normalized) frequency distribution 
P{Es) of the energy of the subsystem, see Fig. [4] 

The hypothesis that Newtonian dynamics causes the sub- 
system to visit points in phase-space with frequencies that 
match with those of the canonical probability distribution can 
now be tested as follows. According to statistical mechanics, 
the distribution of energy in the canonical ensemble is given 
byH 



p(E)=g(E)e- E l T /Z 



(4) 



where T, E, g(E) and Z are the temperature (in units of 
kg = 1), the total energy, the density of states and the partition 
function, respectively |2fl. The function p(E) has a maximum 
at some energy E* , the most probable energy at the tempera- 
ture T [2[|. In the vicinity of E* , we have yfl 



p(E) =Ae' 



a 2 (E-E*) 2 +a 3 {E-E*) 3 +a 4 (E-E*) 3 + .. 



(5) 



where A is a normalization constant and the coefficients a n = 
(1 /n\)d"S(E) I dE"\E* are determined by the microcanonical 
entropy S(E) of the subsystem J2]. 

For the two models considered in this paper, the coefficients 
a n are simple functions of T and Ns (see Supplementary In- 
formation). Therefore, using T and E* as adjustable parame- 
ters a fit of Eq. (O to the histogram P(E S ) obtained from the 
dynamical evolution of the subsystem yields an estimate of 
the temperature Ts of the subsystem. As shown in Fig.|H the 
simulation data for P{Es) (red lines) and fitted p(E) (black 
lines) are in excellent agreement, for both models alike. In 
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FIG. 4. The distribution P{E$) as a function of the energy per par- 
ticle Es/Ns for different sizes Ns of the subsystem, as obtained from 
the solution of Newtonian equations of motion for the oscillator sys- 
tem. Red lines: N s = 20,40,60,80, 100, 120, 160,200 (broad to nar- 
row) and N = 65536. Black lines: probability distribution predicted 
by statistical mechanics Eq. l|5} using all terms up to {E — £*) 8 . In all 
cases, the initial energy of the environment corresponds to a tempera- 
ture T = 1 and the number of samples is 4 x 10 6 . The inset shows the 
corresponding data for the spin system for N$ = 20,40,60,80, 100 
and N = 2000, with the initial energy of the environment correspond- 
ing to T = 0.329 and the number of samples is 10 6 . 

the thermodynamic limit (Ne — > °° before Ns — > °°), all but 
the quadratic term in the exponential can be neglected and the 
distribution is Gaussian Therefore, for large Ns, T$ and 
E* can be obtained by fitting a Gaussian to P(Es) but from 
Fig. 2] it is clear that for small subsystem sizes, P(E S ) de- 
viates significantly from a Gaussian. However, taking into 
account the higher-order terms in the expansion Eq. (0, the 
agreement between simulation data and the prediction of sta- 
tistical mechanics is excellent. Repeating the simulations with 
different initial conditions (including different initial energies 
for the subsystem or the environment) strongly suggests that 
this agreement is generic. 

If the estimate T$ is indeed the temperature of the subsys- 
tem, the second central moment of P{Es) should be related to 
the specific heat of the subsystem. To check this, we define 



where (Es) and are the first and second moment of 

P(Es), respectively. Our numerical results (see Supplemen- 
tary Information) are in excellent agreement with canonical 
ensemble theory. 

As a conclusive test, we perform microcanonical Monte 
Carlo simulations for both models (see Supplementary Infor- 
mation) and obtain the distributions P mc (E s ) of the subsys- 
tems. The microcanonical Monte Carlo simulation generates 
statistically independent configurations of the whole system 
strictly according to the microcanonical distribution but sam- 
ples the phase space in a completely different manner than 
does Newtonian dynamics. The Kullback-Leibler distance 



D(P mc (E s );P(E s )) is a convenient measure to quantify the dif- 
ference between the two distributions P mc {E s ) and P(E S ) l25ll . 
In all cases, we find that the difference between P mc (E s ) and 
P(E S ) is very small (see Supplementary Information). For in- 
stance, for the system of oscillators with Ns = 20, N = 65536 
and 10 8 samples, we find that D(P mc (E s );P(E s )) w 4 x 1(T 2 , 
indicating that the probability that the two distributions are 
different is very small. 

Having established that the interaction of the environment 
with the subsystem causes both systems to equilibrate and 
also drives the latter to its canonical state, it becomes pos- 
sible to derive from the Newtonian dynamics alone, estimates 
for the equilibration time. To this end we express the equili- 
bration time estimated from the simulations in physical units. 
Typical frequencies of vibration in a solid are of the order of 
10 Hz. Using this number to set the scale of the frequency 
Q. in our model, we find that equilibration takes of the order 
of 10~ 9 s. Similarly, for the system of magnetic moments, a 
realistic value of J /kg is of the order of 10 K, yielding an equi- 
libration time of the order of 10~ 8 s. Classical spin systems 
with Hamiltonians that encode frustration and/or disorder of 
regular or random kind are however expected to exhibit larger, 
possibly much larger equilibration time scales. The dynamical 
properties for subsystems of such theories under Newtonian 
evolution are beyond the scope of the present work. 

Even though the subsystems and the environments which 
we have simulated are very small in the thermodynamic sense, 
the subsystem and environment equilibrate on a nanosecond 
time scale. Therefore, for an isolated nanoparticle of even a 
few thousand atoms, an experimental probe that concentrates 
on only a few of those atoms should yield data that follows 
the canonical distribution. 
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